Introduction to Ambisonics¶

Fundamentals¶

Chris Hold, Virtual Acoustics (Aalto University)

Structure¶

  • Motivation
  • Discovering Spherical Harmonics
  • Ambisonics Encoder and Decoder

From Recording to Playback¶

Setup2Playback

Audio Description Formats¶

Channel Based Object Based Scene Based
Channel Object Scene

Scene Based Workflow¶

Channel

Benefits of Scene-based¶

  • Separation of recording, transmission / storage, and playback
  • Flexible with multiple options for each stage
  • Does not scale with the number of sources

Ambisonics¶

  • Implementation of a scene based format
  • "Long" history and connection to other sciences
  • Use spherical harmonics as basis functions

Example¶

No description has been provided for this image
Error in callback <function flush_figures at 0x15c95d1c0> (for post_execute), with arguments args (),kwargs {}:
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib_inline/backend_inline.py:126, in flush_figures()
    123 if InlineBackend.instance().close_figures:
    124     # ignore the tracking, just draw and close all figures
    125     try:
--> 126         return show(True)
    127     except Exception as e:
    128         # safely show traceback if in IPython, else raise
    129         ip = get_ipython()

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib_inline/backend_inline.py:90, in show(close, block)
     88 try:
     89     for figure_manager in Gcf.get_all_fig_managers():
---> 90         display(
     91             figure_manager.canvas.figure,
     92             metadata=_fetch_figure_metadata(figure_manager.canvas.figure)
     93         )
     94 finally:
     95     show._to_draw = []

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/IPython/core/display_functions.py:298, in display(include, exclude, metadata, transient, display_id, raw, clear, *objs, **kwargs)
    296     publish_display_data(data=obj, metadata=metadata, **kwargs)
    297 else:
--> 298     format_dict, md_dict = format(obj, include=include, exclude=exclude)
    299     if not format_dict:
    300         # nothing to display (e.g. _ipython_display_ took over)
    301         continue

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/IPython/core/formatters.py:179, in DisplayFormatter.format(self, obj, include, exclude)
    177 md = None
    178 try:
--> 179     data = formatter(obj)
    180 except:
    181     # FIXME: log the exception
    182     raise

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/decorator.py:232, in decorate.<locals>.fun(*args, **kw)
    230 if not kwsyntax:
    231     args, kw = fix(args, kw, sig)
--> 232 return caller(func, *(extras + args), **kw)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/IPython/core/formatters.py:223, in catch_format_error(method, self, *args, **kwargs)
    221 """show traceback on failed format call"""
    222 try:
--> 223     r = method(self, *args, **kwargs)
    224 except NotImplementedError:
    225     # don't warn on NotImplementedErrors
    226     return self._check_return(None, args[0])

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/IPython/core/formatters.py:340, in BaseFormatter.__call__(self, obj)
    338     pass
    339 else:
--> 340     return printer(obj)
    341 # Finally look for special method names
    342 method = get_real_method(obj, self.print_method)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/IPython/core/pylabtools.py:152, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    149     from matplotlib.backend_bases import FigureCanvasBase
    150     FigureCanvasBase(fig)
--> 152 fig.canvas.print_figure(bytes_io, **kw)
    153 data = bytes_io.getvalue()
    154 if fmt == 'svg':

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/backend_bases.py:2187, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2183 try:
   2184     # _get_renderer may change the figure dpi (as vector formats
   2185     # force the figure dpi to 72), so we need to set it again here.
   2186     with cbook._setattr_cm(self.figure, dpi=dpi):
-> 2187         result = print_method(
   2188             filename,
   2189             facecolor=facecolor,
   2190             edgecolor=edgecolor,
   2191             orientation=orientation,
   2192             bbox_inches_restore=_bbox_inches_restore,
   2193             **kwargs)
   2194 finally:
   2195     if bbox_inches and restore_bbox:

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/backend_bases.py:2043, in FigureCanvasBase._switch_canvas_and_return_print_method.<locals>.<lambda>(*args, **kwargs)
   2039     optional_kws = {  # Passed by print_figure for other renderers.
   2040         "dpi", "facecolor", "edgecolor", "orientation",
   2041         "bbox_inches_restore"}
   2042     skip = optional_kws - {*inspect.signature(meth).parameters}
-> 2043     print_method = functools.wraps(meth)(lambda *args, **kwargs: meth(
   2044         *args, **{k: v for k, v in kwargs.items() if k not in skip}))
   2045 else:  # Let third-parties do as they see fit.
   2046     print_method = meth

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/backends/backend_svg.py:1339, in FigureCanvasSVG.print_svg(self, filename, bbox_inches_restore, metadata)
   1334 w, h = width * 72, height * 72
   1335 renderer = MixedModeRenderer(
   1336     self.figure, width, height, dpi,
   1337     RendererSVG(w, h, fh, image_dpi=dpi, metadata=metadata),
   1338     bbox_inches_restore=bbox_inches_restore)
-> 1339 self.figure.draw(renderer)
   1340 renderer.finalize()

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/artist.py:95, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     93 @wraps(draw)
     94 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 95     result = draw(artist, renderer, *args, **kwargs)
     96     if renderer._rasterizing:
     97         renderer.stop_rasterizing()

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     69     if artist.get_agg_filter() is not None:
     70         renderer.start_filter()
---> 72     return draw(artist, renderer)
     73 finally:
     74     if artist.get_agg_filter() is not None:

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/figure.py:3154, in Figure.draw(self, renderer)
   3151         # ValueError can occur when resizing a window.
   3153 self.patch.draw(renderer)
-> 3154 mimage._draw_list_compositing_images(
   3155     renderer, self, artists, self.suppressComposite)
   3157 for sfig in self.subfigs:
   3158     sfig.draw(renderer)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/image.py:132, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    130 if not_composite or not has_images:
    131     for a in artists:
--> 132         a.draw(renderer)
    133 else:
    134     # Composite any adjacent images together
    135     image_group = []

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/artist.py:39, in _prevent_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     36     renderer.stop_rasterizing()
     37     renderer._rasterizing = False
---> 39 return draw(artist, renderer, *args, **kwargs)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/projections/polar.py:1036, in PolarAxes.draw(self, renderer)
   1033     self.yaxis.reset_ticks()
   1034     self.yaxis.set_clip_path(self.patch)
-> 1036 super().draw(renderer)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     69     if artist.get_agg_filter() is not None:
     70         renderer.start_filter()
---> 72     return draw(artist, renderer)
     73 finally:
     74     if artist.get_agg_filter() is not None:

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/axes/_base.py:3034, in _AxesBase.draw(self, renderer)
   3031     for spine in self.spines.values():
   3032         artists.remove(spine)
-> 3034 self._update_title_position(renderer)
   3036 if not self.axison:
   3037     for _axis in self._axis_map.values():

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/axes/_base.py:2978, in _AxesBase._update_title_position(self, renderer)
   2976 top = max(top, bb.ymax)
   2977 if title.get_text():
-> 2978     ax.yaxis.get_tightbbox(renderer)  # update offsetText
   2979     if ax.yaxis.offsetText.get_text():
   2980         bb = ax.yaxis.offsetText.get_tightbbox(renderer)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/axis.py:1336, in Axis.get_tightbbox(self, renderer, for_layout_only)
   1333     renderer = self.figure._get_renderer()
   1334 ticks_to_draw = self._update_ticks()
-> 1336 self._update_label_position(renderer)
   1338 # go back to just this axis's tick labels
   1339 tlb1, tlb2 = self._get_ticklabel_bboxes(ticks_to_draw, renderer)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/axis.py:2609, in YAxis._update_label_position(self, renderer)
   2605     return
   2607 # get bounding boxes for this axis and any siblings
   2608 # that have been set by `fig.align_ylabels()`
-> 2609 bboxes, bboxes2 = self._get_tick_boxes_siblings(renderer=renderer)
   2610 x, y = self.label.get_position()
   2611 if self.label_position == 'left':

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/axis.py:2161, in Axis._get_tick_boxes_siblings(self, renderer)
   2159 axis = ax._axis_map[name]
   2160 ticks_to_draw = axis._update_ticks()
-> 2161 tlb, tlb2 = axis._get_ticklabel_bboxes(ticks_to_draw, renderer)
   2162 bboxes.extend(tlb)
   2163 bboxes2.extend(tlb2)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/axis.py:1315, in Axis._get_ticklabel_bboxes(self, ticks, renderer)
   1313 if renderer is None:
   1314     renderer = self.figure._get_renderer()
-> 1315 return ([tick.label1.get_window_extent(renderer)
   1316          for tick in ticks if tick.label1.get_visible()],
   1317         [tick.label2.get_window_extent(renderer)
   1318          for tick in ticks if tick.label2.get_visible()])

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/axis.py:1315, in <listcomp>(.0)
   1313 if renderer is None:
   1314     renderer = self.figure._get_renderer()
-> 1315 return ([tick.label1.get_window_extent(renderer)
   1316          for tick in ticks if tick.label1.get_visible()],
   1317         [tick.label2.get_window_extent(renderer)
   1318          for tick in ticks if tick.label2.get_visible()])

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/text.py:958, in Text.get_window_extent(self, renderer, dpi)
    956 bbox, info, descent = self._get_layout(self._renderer)
    957 x, y = self.get_unitless_position()
--> 958 x, y = self.get_transform().transform((x, y))
    959 bbox = bbox.translated(x, y)
    960 return bbox

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/transforms.py:1495, in Transform.transform(self, values)
   1492 values = values.reshape((-1, self.input_dims))
   1494 # Transform the values
-> 1495 res = self.transform_affine(self.transform_non_affine(values))
   1497 # Convert the result back to the shape of the input values.
   1498 if ndim == 0:

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/_api/deprecation.py:297, in rename_parameter.<locals>.wrapper(*args, **kwargs)
    292     warn_deprecated(
    293         since, message=f"The {old!r} parameter of {func.__name__}() "
    294         f"has been renamed {new!r} since Matplotlib {since}; support "
    295         f"for the old name will be dropped %(removal)s.")
    296     kwargs[new] = kwargs.pop(old)
--> 297 return func(*args, **kwargs)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/transforms.py:2419, in CompositeGenericTransform.transform_non_affine(self, values)
   2417     return self._a.transform_non_affine(values)
   2418 else:
-> 2419     return self._b.transform_non_affine(self._a.transform(values))

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/_api/deprecation.py:297, in rename_parameter.<locals>.wrapper(*args, **kwargs)
    292     warn_deprecated(
    293         since, message=f"The {old!r} parameter of {func.__name__}() "
    294         f"has been renamed {new!r} since Matplotlib {since}; support "
    295         f"for the old name will be dropped %(removal)s.")
    296     kwargs[new] = kwargs.pop(old)
--> 297 return func(*args, **kwargs)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/transforms.py:2417, in CompositeGenericTransform.transform_non_affine(self, values)
   2415     return values
   2416 elif not self._a.is_affine and self._b.is_affine:
-> 2417     return self._a.transform_non_affine(values)
   2418 else:
   2419     return self._b.transform_non_affine(self._a.transform(values))

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/_api/deprecation.py:297, in rename_parameter.<locals>.wrapper(*args, **kwargs)
    292     warn_deprecated(
    293         since, message=f"The {old!r} parameter of {func.__name__}() "
    294         f"has been renamed {new!r} since Matplotlib {since}; support "
    295         f"for the old name will be dropped %(removal)s.")
    296     kwargs[new] = kwargs.pop(old)
--> 297 return func(*args, **kwargs)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/transforms.py:2419, in CompositeGenericTransform.transform_non_affine(self, values)
   2417     return self._a.transform_non_affine(values)
   2418 else:
-> 2419     return self._b.transform_non_affine(self._a.transform(values))

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/_api/deprecation.py:297, in rename_parameter.<locals>.wrapper(*args, **kwargs)
    292     warn_deprecated(
    293         since, message=f"The {old!r} parameter of {func.__name__}() "
    294         f"has been renamed {new!r} since Matplotlib {since}; support "
    295         f"for the old name will be dropped %(removal)s.")
    296     kwargs[new] = kwargs.pop(old)
--> 297 return func(*args, **kwargs)

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/matplotlib/projections/polar.py:78, in PolarTransform.transform_non_affine(self, values)
     76     r = (r - self._get_rorigin()) * self._axis.get_rsign()
     77 r = np.where(r >= 0, r, np.nan)
---> 78 return np.column_stack([r * np.cos(theta), r * np.sin(theta)])

File ~/opt/anaconda3/envs/slides/lib/python3.11/site-packages/numpy/lib/shape_base.py:652, in column_stack(tup)
    650         arr = array(arr, copy=False, subok=True, ndmin=2).T
    651     arrays.append(arr)
--> 652 return _nx.concatenate(arrays, 1)

KeyboardInterrupt: 
No description has been provided for this image
No description has been provided for this image

From the Wave Equation to Spherical Harmonics¶

Wave equation \begin{equation} \nabla^2 p - \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2} = 0 \quad. \end{equation}

In frequency domain, with wave-number $k = \frac{\omega}{c}$ results in the Helmholtz equation \begin{equation} (\nabla^2 + k^2) p = 0 \quad. \end{equation}

Multiple solutions, e.g. a mono chromatic plane wave with amplitude $\hat{A}(\omega)$ in Cartesian coordinates \begin{equation} p(t, x, y, z) = \hat{A}(\omega) e^{i(k_x x + k_y y + k_z z - \omega t)} \quad. \end{equation}

No description has been provided for this image

Any solution to the Helmholtz equation can also be expressed in spherical coordinates \begin{equation} p(r, \theta, \phi, k) = \sum_{n = 0}^{\infty} \sum_{m=-n}^{n} \color{olive}{(A_{mn} j_n(kr) + B_{mn} y_n(kr))} \; \color{brown}{Y_{n}^{m}(\theta,\phi)} \quad, \end{equation}

\begin{equation} p(r, \theta, \phi, k) = \sum_{n = 0}^{\infty} \sum_{m=-n}^{n} \color{olive}{C_{mn}(kr)} \; \color{brown}{Y_{n}^{m}(\theta,\phi)} \quad, \end{equation}

With two separable parts:

  • radial component : linear combination of spherical Bessel functions of first ($j_n$) and second kind $y_n$
  • angular component : spherical harmonics $\color{brown}{Y_{n}^{m}(\theta,\phi)}$

sound field defined on a sphere is fully captured by its spherical harmonics coefficients

E.g. a unit plane wave \begin{equation} p(r, \theta, \phi, k) = \sum_{n = 0}^{\infty} \sum_{m=-n}^{n} 4 \pi i^n j_n(kr) \left[ Y_n^m(\theta_k,\phi_k) \right] ^* Y_{n}^ {m}(\theta,\phi) \quad. \end{equation}

Discovering Spherical Harmonics¶

\begin{equation} Y_n^m(\theta,\phi)=\color{orange}{\sqrt{\frac{2n+1}{4\pi}\frac{(n-m)!}{(n+m)!}}}\, \color{teal}{P_n^m(\cos(\theta))}\, \color{purple}{e^{im\phi}} \quad, \end{equation}\begin{equation} Y_n^m(\theta,\phi)=\color{orange}{D_{nm}}\, \color{teal}{P_n^m(\cos(\theta))}\, \color{purple}{e^{im\phi}} \quad, \end{equation}

where

  • azimuth as $\phi$ and zenith/colatitude as $\theta$.
  • $P_n^m$ is the associated Legendre polynomial of order $n$ and degree $m$.

Combined with appropriate scaling give real Spherical Harmonics $Y_{n,m}(\theta,\phi)$ as

\begin{equation} Y_{n,m}(\theta,\phi) = \sqrt{ \frac{(2n+1)}{4\pi} \frac{(n-|m|)!}{(n+|m|)!} } P_{n,|m|}(\cos(\theta)) \begin{cases} \sqrt2\sin(|m|\phi) & \mathrm{if\hspace{0.5em}} m < 0 \quad,\\ 1 & \mathrm{if\hspace{0.5em}} m = 0 \quad,\\ \sqrt2\cos(|m|\phi) & \mathrm{if\hspace{0.5em}} m > 0 \quad. \end{cases} \end{equation}
  • Always check: Condon–Shortley phase convention $(-1)^m$
Error in callback <function _draw_all_if_interactive at 0x11220e520> (for post_execute), with arguments args (),kwargs {}:
KeyboardInterrupt

No description has been provided for this image

Let's look at the azimuthal component $e^{im\phi}$ and the zenithal component $P_n^m(\cos\theta)$

azimuthal component¶

No description has been provided for this image
No description has been provided for this image

zenithal component¶

No description has been provided for this image

Orthonormality¶

Two functions $\color{violet}f, \color{purple}g$ over a domain $\gamma$ are orthogonal if

\begin{equation} \int_\gamma \color{violet}{f^*(\gamma)} \color{purple}{g(\gamma)} \,\mathrm{d}\gamma = \langle \color{violet}f, \color{purple}g \rangle = 0 \quad,\,\mathrm{for}~f \neq g. \end{equation}

They are also orthonormal if \begin{equation} \color{violet}{\int_\gamma f^*(\gamma) f(\gamma) \,\mathrm{d}\gamma = \int_\gamma|f(\gamma)|^2 \,\mathrm{d}\gamma = \langle f, f \rangle } = 1 \quad, \\ \color{purple}{\int_\gamma g^*(\gamma) g(\gamma) \,\mathrm{d}\gamma = \int_\gamma |g(\gamma)|^2 \,\mathrm{d}\gamma = \langle g, g \rangle} = 1 \quad. \end{equation}

For the Spherical Harmonics:

  • the the azimuthal component $e^{im\phi}$ along $\phi$ is orthogonal w.r.t. the degree $m$
  • the zenithal component $P_n^m(\cos\theta)$ along $\theta$ is orthogonal w.r.t. the order $n$.

Their product is still orthogonal, and the scaling $\color{orange}{D_{nm}}$ ensures orthonormality such that \begin{equation} \int_\Omega Y_n^m(\Omega)^* \, Y_{n'}^{m'}(\Omega) \,\mathrm{d}\Omega = \langle Y_n^m(\Omega) , Y_{n'}^{m'}(\Omega) \rangle = \delta_{nn'}\delta_{mm'} \quad , \end{equation}

and

\begin{equation} \int_{{\Omega} \in \mathbb{S}^2} |Y_n^m({\Omega})|^2 \mathrm{d}{\Omega} = 1 \quad . \end{equation}

Example¶

show normality for $Y_0^0$: $$ Y_0^0(\theta,\phi)=\sqrt{\frac{0+1}{4\pi}\frac{(0)!}{(0)!}} P_0^0(\cos(\theta)) e^{i0\phi} = \sqrt{\frac{1}{4\pi}} \quad, $$ hence $$ \int_{{\Omega} \in \mathbb{S}^2} Y_0^0(\theta,\phi)^* Y_0^0(\theta,\phi) \,\mathrm{d}{\Omega} = \int_{{\Omega} \in \mathbb{S}^2} \sqrt{\frac{1}{4\pi}}\sqrt{\frac{1}{4\pi}} \, \mathrm{d}{\Omega} = 4\pi \frac{1}{4\pi} = 1 \quad. $$

Example¶

Test orthogonality of functions $\color{violet}f, \color{purple}g$

$$ \langle \color{violet}f, \color{purple}g \rangle \overset{?}{=} 0 \quad,\,\mathrm{for}~f \neq g. $$
Test for orthogonality in azimuth <f_azi, g_azi> =  -0.0
Test for orthogonality in zenith <f_zen, g_zen> =  0.0
No description has been provided for this image

Spherical Harmonic Transform (SHT)¶

  • Find a scene based encoding for sound field
  • We showed that the spherical harmonics $Y_n^m({\Omega})$ are a set of suitable basis functions on the sphere
  • We also showed that a sound field (on the sphere) $s({\Omega})$ is fully captured by its spherical harmonics coefficients $\sigma_{nm}$

This can be expressed with $\Omega = [\phi, \theta]$ as the inverse Spherical Harmonic Transform (iSHT) \begin{equation} s({\Omega}) = \sum_{n = 0}^{N=\infty} \sum_{m=-n}^{+n} \sigma_{nm} Y_n^m({\Omega}) \quad. \end{equation}

Spherical harmonics coefficients $\sigma_{nm}$ can be derived with the Spherical Harmonic Transform (SHT) \begin{equation} \sigma_{nm} = \int_{{\Omega} \in \mathbb{S}^2} s({\Omega}) [Y_n^m({\Omega})]^* \mathrm{d}{\Omega} = \langle [Y_n^m({\Omega})] , s({\Omega}) \rangle \quad. \end{equation}

  • SHT also refered to as spherical Fourier transform
  • Note the transform from a spatially discrete to continuous domain

Spherical Grids¶

The SHT evaluates the continuous integral over $\Omega$ \begin{equation} \sigma_{nm} = \int_{{\Omega} \in \mathbb{S}^2} s({\Omega}) [Y_n^m({\Omega})]^* \mathrm{d}{\Omega} \quad . \end{equation} Quadrature methods allow evaluation by spherical sampling at certain (weighted) grid points such that \begin{equation} \sigma_{nm} \approx \sum_{q=1}^{Q} w_q s({\Omega}_q) [Y_n^m({\Omega}_q)]^* \quad. \end{equation}

Certain grids with sampling points ${\Omega}_q$ and associated sampling weights $w_q$ have certain properties:

  • quadrature grids allow numerical integration of spherical polynomials
  • Spherical sampling dictates maximum order
  • easy to evaluate are uniform/regular grids, with constant $w_q = \frac{4\pi}{Q}$
  • An example are so-called spherical t-designs(t), which allow a SHT up to order $N = \lfloor t/2 \rfloor$
No description has been provided for this image
No description has been provided for this image

Spherical Dirac¶

We can show that spherical harmonics are orthogonal (even orthonormal) with \begin{equation} \int_\Omega Y_n^m(\Omega) \, Y_{n'}^{m'}(\Omega) \,\mathrm{d}\Omega = \delta_{nn'}\delta_{mm'} \quad . \end{equation}

Because of their completeness, we can also directly formulate a Dirac function on the sphere as \begin{equation} \sum_{n=0}^{N=\infty} \sum_{m=-n}^n [Y_n^m({\Omega'})]^* Y_n^m(\Omega) = \delta(\Omega - \Omega') \quad, \end{equation}

and therefore the spherical Fourier coefficients $\sigma_{nm}$ \begin{equation} SHT\{\delta(\Omega - \Omega')\} = \int_{{\Omega} \in \mathbb{S}^2} \delta(\Omega - \Omega') \, [Y_n^m({\Omega})]^* \mathrm{d}{\Omega} = [Y_n^m({\Omega'})]^* \quad . \end{equation}

Order-Limitation of Spherical Dirac Pulse¶

No description has been provided for this image

Example¶

Integrate (order-limited) Dirac $\delta_N(\Omega - \Omega')$ over sphere

$$ \int_{{\Omega} \in \mathbb{S}^2} \delta_N(\Omega - \Omega') \mathrm{d}{\Omega} = \int_{{\Omega} \in \mathbb{S}^2} \color{blue}{\sum_{n=0}^{N} \sum_{m=-n}^n\, [Y_n^m({\Omega'})]^* Y_n^m(\Omega)} \,\mathrm{d}{\Omega} \quad, $$

by discretization with sufficient t-design $$ \int_{{\Omega} \in \mathbb{S}^2} \color{blue}{\sum_{n=0}^{N} \sum_{m=-n}^n\, [Y_n^m({\Omega'})]^* Y_n^m(\Omega)} \,\mathrm{d}{\Omega} = \frac{4\pi}{Q}\sum_{q=1}^{Q} \color{orange}{\sum_{n=0}^{N} \sum_{m=-n}^n\, [Y_n^m({\Omega'})]^* Y_n^m(\Omega_q)} $$

No description has been provided for this image

Matrix Notations¶

Stack the spherical harmonics evaluated at $\Omega$ up to spherical order $N$ as

$$ \mathbf{Y} = \left[ \begin{array}{ccccc} Y_0^0(\Omega[0]) & Y_1^{-1}(\Omega[0]) & Y_1^0(\Omega[0]) & \dots & Y_N^N(\Omega[0]) \\ Y_0^0(\Omega[1]) & Y_1^{-1}(\Omega[1]) & Y_1^0(\Omega[1]) & \dots & Y_N^N(\Omega[1]) \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ Y_0^0(\Omega[Q-1]) & Y_1^{-1}(\Omega[Q-1]) & Y_1^0(\Omega[Q-1]) & \dots & Y_N^N(\Omega[Q-1]) \end{array} \right] $$

such that $\mathbf{Y} : [Q, (N+1)^2]$

  • soundfield pressure is real signal -> only real spherical harmonics needed
  • orthonormal scaling introduced earlier called N3D by convention (in contrast to SN3D)
  • stacking them as above, where $idx_{n,m} = n^2+n+m$, is called ACN
  • We can stack a time signal $s(t_0, t_1, \ldots, t_T)$ as vector $\mathbf{s}$.
\begin{equation*} \mathbf{s} = \begin{bmatrix} s(t_0) \\ s(t_1) \\ \vdots \\ s(t_T) \end{bmatrix} \end{equation*}
  • This leads to a matrix notation of multiple discrete signals $\mathbf{S} : [T, Q]$.
\begin{equation*} \mathbf{S} = \begin{bmatrix} s_0(t_0) & s_1(t_0) & \cdots & s_Q(t_0)\\ s_0(t_1) & s_1(t_1) & \cdots & s_Q(t_1)\\ \vdots & \vdots& \ddots & \vdots \\ s_0(t_T) & s_1(t_T) & \cdots & s_Q(t_T) \end{bmatrix} \end{equation*}
  • Similarly, stacking Ambisonics coefficients $\sigma_n^m(t_0, t_1, \ldots, t_T)$ into $\mathbf{\chi} : [T, (N+1)^2]$
\begin{equation*} \mathbf{\chi} = \begin{bmatrix} \sigma_0^0(t_0) & \sigma_1^{-1}(t_0) & \cdots & \sigma_N^N(t_0) \\ \sigma_0^0(t_1) & \sigma_1^{-1}(t_1) & \cdots & \sigma_N^N(t_1) \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_0^0(t_T) & \sigma_1^{-1}(t_T) & \cdots & \sigma_N^N(t_T) \end{bmatrix} \end{equation*}

We obtain ambisonic signals matrix $\mathbf{\chi} : [T, (N+1)^2]$ from signals $\mathbf{S} : [T, Q]$ by SHT in matrix notation as \begin{equation} \mathbf{\chi}(t) = \mathbf{S}(t) \, \mathrm{diag}(w_q) \, \mathbf{Y} \quad. \end{equation} With the SH basis functions evaluated at $\Omega_q$ as $\mathbf{Y} : [Q, (N+1)^2]$.

Obtaining the discrete signals $s_q(t)$ is a linear combination of SH basis functions evaluated at $\Omega_q$ as $\mathbf{Y} : [Q, (N+1)^2]$ and the SH coefficients. This inverse SHT in matrix notation with ambisonic signals matrix $\mathbf{\chi} : [T, (N+1)^2]$ is

\begin{equation} s_q(t) = \mathbf{\chi}(t) \, \mathbf{Y}^T \quad . \end{equation}

From the Ambisonics Encoder, to Directional Weighting, to a Decoder

Encoder¶

A single plane wave encoded in direction $\Omega_1$ with signal $\mathbf{s}$ is directly the outer product with the Dirac SH coefficients

\begin{equation} \mathbf{\chi}_\mathrm{PW(\Omega_1)} = \mathbf{s} \, \mathbf{Y}(\Omega_1) \quad. \end{equation}

For multiple sources $Q$, we stack and sum \begin{equation} \mathbf{\chi}_\mathrm{PW(\Omega_Q)} = \sum_{q=1}^Q\mathbf{s}_q \, \mathbf{Y}(\Omega_q) = \mathbf{S} \, \mathbf{Y} \quad. \end{equation}

No description has been provided for this image
No description has been provided for this image

Directional Weighting¶

  • Directionally filtering a soundfield in direction $\Omega_k$ by weighting
  • Weighting in SH domain $w_{nm}$ is an elegant way of beamforming

The simplest beamformer is a spherical Dirac in direction $\Omega_k$ normalized by its energy, i.e. $\mathrm{max}DI$ \begin{equation} w_{nm, \mathrm{maxDI}}(\Omega_k) = \frac{4\pi}{(N+1)^2} Y_{n,m}(\Omega_k) \end{equation}

No description has been provided for this image

Other patterns can be achieved by weighting the spherical Fourier spectrum. Axis-symmetric patterns reduce to only a modal weighting $c_n$, such that \begin{equation} w_{nm}(\Omega_k) = c_{n} \, Y_{n,m}(\Omega_k) \end{equation}

E.g. $\mathrm{max}\vec{r}_E$ weights each order $n$ with \begin{equation} c_{n,\,\mathrm{max}\vec{r}_E} = P_n[\cos(\frac{137.9^\circ}{N+1.51})] \quad, \end{equation} with the Legendre polynomials $P_n$ of order $n$.

No description has been provided for this image

Or we can define a (SH) Butterworth filter with \begin{equation} c_{n,\,\mathrm{Butterworth}} = \frac{1}{\sqrt{1+(n/n_c)^{2k}}} \quad, \end{equation} with the filter-order $k$ and cut-on SH-order $n_c$.

No description has been provided for this image

Decoder¶

  • Beamformer extracts a portion of the soundfield in steering direction $\Omega$
  • Pointing a set of beamformers in directions $\Omega_k$ results in a set of signal components
  • in case of $\mathrm{max}DI$ proportional to iSHT in $\Omega_k$
\begin{equation} s({\Omega_k}) = \sum_{n = 0}^{N} \sum_{m=-n}^{+n} w_{nm}({\Omega_k}) \, \sigma_{nm} \quad, \end{equation}

or in matrix notation with $\mathbf{S}$ and beamforming weights $\mathbf{c}_n$ and $\mathbf{Y}$ \begin{equation} \mathbf{S} = \mathbf{\chi} \, \mathrm{diag_N}(\mathbf{c}_n) \, \mathbf{Y}^T \quad . \end{equation}

Example¶

SHD signal of 3 sine signal PWs plus noise -> three beamformers, two pointing to sources

No description has been provided for this image
No description has been provided for this image

Example¶

Decode on a t-design(6) (sufficient up to $N = 3$):

  • a 3rd order ambisonic signal
No description has been provided for this image

Example¶

Decode on a t-design(6) (sufficient up to $N = 3$):

  • a 5th order ambisonic signal
  • an 8th order ambisonic signal
No description has been provided for this image
No description has been provided for this image

Loudspeaker Decoders¶

In [25]:
# Loudspeaker Setup
ls_dirs = np.array([[-135, -80, -45, 0, 45, 80, 135, -135, -60, -30, 30, 60, 135],
                    [0, 0, 0, 0, 0, 0, 0, 60, 60, 60, 60, 60, 60]])
ls_x, ls_y, ls_z = spa.utils.sph2cart(spa.utils.deg2rad(ls_dirs[0, :]),
                                      spa.utils.deg2rad(90 - ls_dirs[1, :]))
ls_setup = spa.decoder.LoudspeakerSetup(ls_x, ls_y, ls_z)
ls_setup.pop_triangles(normal_limit=85, aperture_limit=90,
                       opening_limit=150)
ls_setup.show()
Face not pointing towards listener: [0 1 2]
Face not pointing towards listener: [0 6 2]
Face not pointing towards listener: [2 5 6]
Face not pointing towards listener: [2 3 4]
Face not pointing towards listener: [2 5 4]
No description has been provided for this image
In [26]:
spa.plot.decoder_performance(ls_setup, 'NLS')
No description has been provided for this image
In [27]:
spa.plot.decoder_performance(ls_setup, 'VBAP')
No description has been provided for this image
In [28]:
spa.plot.decoder_performance(ls_setup, 'SAD')
No description has been provided for this image
In [29]:
spa.plot.decoder_performance(ls_setup, 'MAD')
No description has been provided for this image
In [30]:
spa.plot.decoder_performance(ls_setup, 'EPAD', N_sph=2)
No description has been provided for this image
In [31]:
ls_setup.ambisonics_setup(update_hull=True)
spa.plot.decoder_performance(ls_setup, 'ALLRAD')
No description has been provided for this image
In [32]:
ls_setup.ambisonics_setup(update_hull=True)
spa.plot.decoder_performance(ls_setup, 'ALLRAD2')
No description has been provided for this image

Binaural Decoding¶

\begin{equation} s^{l,r}(t) = x(t) * h_{\mathrm{HRIR}}^{l,r}({\Omega}, t) \quad, \end{equation}

where $(*)$ denotes the time-domain convolution operation.

Transforming to the time-frequency domain through the time-domain Fourier transform, further assuming plane-wave components $\bar X (\Omega)$, the ear input signals are given as \begin{equation} S^{l,r}(\omega) = \int_{\Omega} \bar X (\Omega, \omega) H^{l,r}(\Omega, \omega) \,\mathrm{d}\Omega \quad. \end{equation}

\begin{equation} S^{l,r}(\omega) = \sum_{n = 0}^{N} \sum_{m = -n}^{+n} \chi_{nm}(\omega) \breve H_{nm}^{l,r}(\omega) \quad. \end{equation}

For one ear (left) this can be interpreted as a frequency dependent ambisonic beamformer \begin{equation} s^l(\omega) = \chi_{nm}(\omega) [\breve H_{nm}^{l}(\omega)]^T \quad . \end{equation}

In [ ]: